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SUMMARY 


The  determination  of  the  equilibrium  composition  of  a 
gaseous  mixture  is  equivalent  to  the  determination  of  i  e 
number  of  moles  of  each  molecular  species  present  that  minimize 
the  total  free  energy  of  the  mixture.  The  convex  function 
for  free  energy  as  given  by  Gibbs  is  minimized  subject  to  the 
constraints  of  the  ina3o— balance  equations  for  the  mixture. 

A  solution  to  the  problem  is  given  using  a  version  of  the 
decomposition  procedure  for  linear  programming.  A  mathematical 
description  of  the  method  of  solution  as  well  as  the  compu¬ 
tational  algorithm  as  programmed  for  the  IBM  704  with  some 
results  are  given. 
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SOLVING  THE  CHEMICAL  EQUILIBRIUM  PROBLEM  USING  THE 
DECOMPOSITION  PRINCIPLE 
George  B.  Dantzig  and  Marvin  Shapiro 

An  earlier  paper  by  one  of  the  authors  on  Convex  Objectives 
[5]  provides  the  theoretical  foundation  for  this  report.  It  gave 
a  proof  of  convergence  for  the  case  when  an  infinite  master 
program  is  generated  by  the  decomposition  process  [2].  In  this 
paper  the  method  Is  applied  to  the  convex  function  of  Gibbs 
Whloh  represents  the  free  energy  of  a  chemical  system  under 
oonstant  temperature  and  pressure.  This  procedure,  which  Is 
similar  In  some  respects  to  that  discussed  In  [4] ,  was  compared 
with  the  "quadratic  fit"  procedure  [3]  further  developed  by 
Seiner  Johnson  and  Herschel  Ranter  and  now  In  wide  use.  This 
procedure  was  extended  by  S.  Johnson  and  one  of  the  authors 
to  the  Multlphased  Chemical  Systems  In  [6]  and  the  code  described 
In  Section  *  covers  this  case  also.  In  each  test  problem  tried 
the  quadratic  fit  appeared  to  be  more  efficient.  This  was  true 
even  after  correction  for  the  use  of  single  precision  and  the 
use  of  Internal  memory  for  the  quadratic  fit  technique  Instead 
of  the  use  of  double  precision  and  tape  for  the  decomposition 
approach.  Empirical  experience  at  this  point  Is  too  limited 
to  say  more  than  that  on  larger  problems  the  two  methods 
seemed  more  oomparable.  For  example  a  problem  Involving  12 
mass  balance  constraints  and  30  molecular  species,  see  [6], 
took  17  Iterations  which  Is  equivalent  to  17  x  12  -  204  Iterations 
of  the  decomposition  principle;  actually  the  latter  required  230. 
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However,  the  present  method  has  certain  advantages.  It  is  more 
flexible  because,  with  minor  code  changes,  it  can  be  extended 
(as  will  be  developed  in  a  future  report)  to  solve  a  very  broad 
class  of  problems;  namely  to  find 

Min  0Q(x)  x  -  (x1,x2,...,xn)  . 

subject  to  m  inequality  constraints  0^(x)  <[  0  where 

rJ1(x)  for  i  =  0,  1,  2,...,m  are  general  convex  functions  in 

n  variables  x  =  (x^,x2> • • • >xn)  over  a  closed  bounded  convex 


domain  x  e  R. 
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1_. _ GENERAL  THEORY 

We  oonsider  an  equilibrium  mixture  containing  m  different 
atom  types.  While  in  theory  these  will  combine  into  all  chemi¬ 
cally  possible  molecular  species,  in  practice  only  standard 
types  are  considered,  including  the  monotonic  types,  which 
are  known  to  occur  in  measurable  amounts. 

Let 

b^  -  the  number  of  atomic  weights  of  species  i 
present  in  the  mixture, 

Xj  -  the  number  of  moles  of  molecular  species  J 
present  in  the  mixture  where 

(1)  xj  2  °'  J  -  1,  2,  ...,  n, 

x  -  the  total  number  of  moles  of  gas  in  the 
mixture,  l.e., 

(2) 

a^  «  the  number  of  atoms  of  species  i  in  a 
molecule  of  species  J. 

Then  the  mass  balance  equations  are 
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The  determination  of  the  equilibrium  composition  of  a 
gaseous  mixture  is  equivalent  to  the  determination  of  the 
values  of  the  mole  numbers  x^  that  obey  constraints  and 
minimize  the  total  free  energy  of  the  mixture  given  by 

n  n 

(4)  G(x1#  . ..,  xn)  -£  °jxj  +  £  xj  log 

J-l  J-l 


n 

L 

j-i 


n 

°j(x/x“  )  +  £  (Xj/x) 
J-l 


log  (x  /x) 


Xj/x  2:  0. 

which  can  be  shown  to  be  a  convex  function.  The  values  0 

J 

are  modified  Qibba  free  energy  funotion  F°/RT  of  the  atomio 
speoies  at  a  given  temperature  plus  the  natural  logarithm  of 
the  pressure  in  atmospheres.  Our  problem  is  to  minimize 

(4)  subject  to  the  linear  equality  and  inequality  constraints 

(1),  (2),  (}). 

Homogeneous  objectives  of  degree  I»  In  the  chemical 
equilibrium  problem  G(x)  appears  in  the  speoial  form 

(5)  0(x)  -  xO(lL,  *4,  .... 

xx  x 

where 

(6)  x  -  x1  +  x2  +  ...  +  xni 

that  is  to  say,  0  is  a  homogeneous  funotion  of  degree  1. 


-  x  £  (o,u,  +  u,  log  u4),  u.  - 


J  J 


is  easily  seen  to  be  a  convex  function. 

In  this  case,  we  may  rewrite  our  system  (2)  in  the  form, 
find  x  ^  0,  and  ratios  u^  -  x^/x  £  0,  Min  z  satisfying 


(7)  x 


a  Uj  +  ai2u2  +  ***  +  a 


lnun 


•  D. 


-  z  ( Min ) 

-  1,  u,  >  0. 

j  L- 

variable  (x)  m—equation 
linear  programming  problem  with  variable  coefficients 

(9)  ^  ^ ‘ ' 1 

y  >  o(u), 

where  u  satisfies  (8).  It  is  easy  to  prove  that  if  0(u) 
is  a  convex  function,  the  set  of  coefficients  [a^,  a0,...,a  ;  y] 
are  selected  from  a  convex  set.  Henoe  we  may  apply  the 
generalized  coefficient  method  discussed  in  The  RAND  Corporation 
Paper  P— 15^4  on  the  Decomposition  Principle  [2]  and  the 


a_nun  +  a_^u^  +  ...  +  a _ u^ 


mil  m2  2 


xO(u) 


mn 


where 

(8) 


u,  +  u0  +  . .  .  +  u, 

jl  c  n 


We  regard  system  (7)  as  a  single 
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special  application  discussed  in  P—1664  on  General  Convex 
Objective  Forma  [J5  J  . 

Iterative  procedure:  We  assume  the  algorithm  can  be 
initiated  with  some  m  sett  j£  veotors  vP  satisfying  (8), 
which  generate 

(10)  alp  -  £  aAJup  p  -  1,  2,  m 

J 

Tp  -  O(uP). 

such  that  Caip]  forms  a  feasible  basis  for  representing  b,  i.e. 
there  exist  *p  -  £  0,  *  ■  Zq  8atisfying 

Ul)  Z  alpXp  "  bi'  p  "  2>  "•>  m 

P 

7  y  \  -  z. 

L  PP 


The  simplex  multipliers 


tt°  satisfy 


(12) 


z 

1 


v.a. 
i  iP 


p  **  1  f  2,  •  •  •  p  m  * 


To  test  optimality,  the  expression 


(13)  0(u)  -  J  ir°  2  aijuj  "  A 

is  minimized  over  all  u^  >  0  satisfying 
Min  A  >  0,  the  solution 


1.  If 


(14)  x-^Upx° 

P 
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r 

rry-l-  1 

is  optimum.  If  not,  choose  u  -  u  such  that 
(15)  0(  um+1)  -  l  l  a1jUjin+1  -  Min  A. 

Augment  system  ( 11 )  and  find  Xp  -  X^  >  0,  Min  z  -  "z ^  satisfying 

^l6)  E  aipxp  +  aim+lVl-l  "  bi 
E  Tpxp  +  VtlNn+l  “  2 

where  ^m+1  are  **ound  by  substituting  u  in 

(10).  The  new  simplex  multipliers  are  substituted  for 
ir°  in  expression  (13)  andA  is  minimized  over  all  £  u^  -  1,  u^  £  0. 
The  procedure  is  then  iterated  introducing  more  and  more  oolumns 
m+1,  m+2,  ...  into  (l6). 

Theorem!  Given  a  homogeneous  convex  objective  form  0(x) 
of  degree  unity  such  that  Inf  0(x)  exists ,  then  the  modified 
algorithm  converges  to  an  optimal  solution  under  e  —perturbation 
if  there  exists  at  least  one  nondegenerate  baslo  solution  to 
the  "master"  program  (l6). 

For  this  application,  then, 

(17)  A  -  o(u)  -  E  2  auuj 
-  Z  (Uj  log  “j  +  Vj) 

A  proof  of  this  theorem  will  be  found  in  The  RAND  Corporation 
P—1664.  It  is  assumed  Inf  G(x)  exists,  where  x  is  in  a  convex  set 
defined  by  (l)  and  (5). 
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where  -  tt^  are  the  simplex  multipliers  of  some  iteration 
and 

m 

(is)  oj  -  oj  - 1  vir 

i-1 

To  find  the  Min  A  subject  to 


(19)  Z  uj  -  i* 


UJ^°' 


we  Ignore  the  relations  u 
minimum  of  the  function 


^  0  and  find  the  unconditional 


(2°)  ^  “  Z  (uj  loS  uj  +  °juj)  ~  6  (Z  uj  “ 

where  0  is  a  Lagrange  multiplier.  We  set  the  partial 
derivatives  of  A  with  respect  to  Uj  equal  to  zero;  thus 

(21)  -  0  -  1  —  0  +  log  u.  +  o’., 

7)Uj  J  J 


whence  u^  may  be  written  in  the  form 
— c  . 

(22)  Uj  -  Ae  J 


where  A  »  e  >  0.  Substituting  into  (19)  determines 
A  >  0  so  that  the  conditions  u^  £  0  hold  at  the  minimum; 

hence 


(2J)  - 


>  0. 


e 
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k  +1 

These  are  the  values  u^  »  to  initiate  the  next  iteration. 

The  value  of  Min  A  is  obtained  by  substituting  (2J)  into 

(i7): 

Min  A  »  log  A  *»  -  log  £  e  d 


Letting  x  >  0  be  any  vector  such  that  Ax  -  b,  then 

Q(x)  _  irb  »  G(x)  —  irAx  -  x  [p(u)  —  ttAuT)  £  X  Min  A 


so  that  Trb  because  of  convexity  and  definition  of  ir  provides 
an  upper  bound  for  G(x)  and  irb  —  x  Min  A  provides  a  lower 
bound  where  x  is  the  optimum  x  and  x  1b  the  sum  of 
components  of  x.  Hence  Min  A  may  be  used  as  a  measure  of 
the  maximum  possible  decrease  in  free  energy  per  mole. 


2.  COMPUTATIONAL  ALGORITHM 


Find  x  >  0  and  Min  G(x)  satisfying 

J 


+  a,  . ,  x  a,  x  -  b. 

lm+1  m+1  In  n  1 


x  +  a  , t  x  ,,+...+  a  x  -b 
m  mm+1  m+1  mn  n  m 


E  Xj  (cj  +  in  Xj/x)  -  0(X) 


where  x  »  x^  +  +  ...  +  xn  and  >  J 
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We  define  A  -  [a^j] ,  b  -  (b^,  bg» • • • ,bm) ,  column  vector; 
x  -  (x1#  Xg ,...,xn),  column  vector; 

Ax  -  o,  column  vector. 

Each  Iteration  p  begins  by  computing  (l)  a  "candidate" 
x«xp;  (2)  Axp-  ap;  (3)  o(xp)  -  y p.  A  test  Is  then  made  to  see 
how  close  Min  A  Is  to  zero.  If  MlnA»Whlch  is  a  measure  of  the 
non-optlmallty  of  the  solution.  Is  within  specified  limits,  say 
h*  <  0,  the  "solution"  Is  printed  out  and  the  computation 
terminates.  The  Information  on  hand  at  the  start  of  the  p^1 
Iteration  consists  of  m+3  columns  of  Information  of  which  m 
are  associated  with  some  subset  k-k1#  k2#***,km  °*  m  PrevlouB 

\r 

x  Including  a  special  Initiating  set.  For  notatlonal  convenience 
below  we  will  assume  k-1,  2,...,m  and  p  )  m  but  will  understand 
that  k  actually  takes  on  some  se*  of  values  k^.  Let 
ak-(olk,...,amk).  Then  the  data  at  hand  at  the  start  of  an 
Iteration  can  be  obtained  from  the  "Information  matrix. " 


1  0  1  n  y2  ...  ym 

0 

0  1  1  1  1  ...  1 
_  | 

0 

0  0  j  «u  a12... 

0  0  |  a21  a22*  °2m 

e 

CM 

1  : 

e 

e 

0  0  |  aml  °m2  *  *  *  amm 

bm 
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The  square  array  B-  (a  Is  called  the  basis  associated  with 
the  iteration.  It  has  a  special  property,  namely,  if  we 

determine  a  column  vector  such  that  BX  -  b.  Its  components 
X^  2.  0.  Instead  of  recording  explicitly  the  above  matrix  the 
set  of  ’'inverse"  values  are  recorded 


where 


i —  These  values  known 
at  3 tart  of  Cycle 


and_ 

x* 

y 


p 


The  above  relations  are  not  used  in  the  computation;  they  are 
shown  to  give  the  relation  between  the  "information"  matrix 
and  its  "inverse." 

To  initiate  the  algorithm  the  special  lnltatlng  set  is 


“  ( 1 ,  0,  0,  ....  0) 
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x02  -  (0,  1,  0,  ....  0) 

x0rn  -  (0 -  1,  0 - 0) 

Thus  the  starting  basis  Is  the  Identity  In  the  "information 
matrix”  below 


1 

o  1 

C, 

. . .  .  c 

0 

0 

1 

1 1 

1 

1 

1  ... 

m 

.  .  .  .  1 

0 

0 

0  1 

1 

0  ... 

.  .  .  0 

bl 

0 

0  1 

0 

1  ... 

...  0 

b2 

0 

1 

1 

o  i 

0 

• 

• 

• 

0  ... 

...  1 

• 

• 

• 

b 

_ 

I 

1 

m 

In  inverse  form, to  start  cycle  1  we  have: 


ITERATIVE  PROCEDURE 
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m 

Step  I:  Compute  c  ,  =  c  .  -  £  a^  1 


— c  . 

yj  ■  e  J: 


y 

y  v  cr°  —Iny 


(j=l,2,...,n) 

y  -  Z 

T 

y  -  (y1>  y2*  • • • >vn) 


Step  II:  Adjoin  a  new  column  to  the  set  of  Inverse  values 


1 

0 

i-7ri 

— TTn  .  .  . 

— 7T 

m 

“g 

h 

0 

1 

j 

|p01 

_^0m_  _ 

— d 

e 

0 

0 

iP , 
l  R 

0 1 2  •  •  • 

Plm 

h 

*1 

0 

0 

iPpi 

1  ; 

Pr2  ‘  * 

Prm 

K 

© 

0 

0 

yi 

Pm2  * • • 

Pim 

where 

h  =  — +  y3 

e  =>  3  a^+1 

0 

M-  c  [6jLj-^a^ 


—iny 


Choose  r  so  that 


where  u^,  p.  )  0  and  r,  1  =  1, 2,..., m 
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Step  III:  If  y  «  1  (within  an  agreed  upon  error  range) 
print  out  solution 

x  »  ^ l  xi  i  D  1 ,  . . . #  m 

t  h 

where  x  is  the  solution  x  which  generated  l  column 
a  i  in  the  basis.  If  y  +  1  go  to  step  IV. 


Step  IV:  Pivot  on  |i^ 

Step  V:  Initiate  new  iteration  with  transformed  tableau 
dropping  adjoined  column. 


3.  ALTERNATIVE  COMPUTATIONAL  PROCEDURE 

The  "information  matrix"  at  the  start  of  the  p  iteration 
is  the  same  as  before  except  that  it  is  arranged  in  the  form 


1 

0 

0 

1 

1 

1 

•  •  •  • 

0 

0 

1 

0 

Y1 

Y  2  •  •  • 

....  y 

m 

0 

0 

0 

1 

1 

1 

•  •  •  •  1 

0 

—  - 

-  —  - 

-  —  J 

_______ 

•  — 

0 

0 

0 

a, b. 

b. 

0 

0 

» 

0 

11  1 

a21"b2 

• 

12  1 

•  •  • 

’ ’ * '  1ml 
-  a2m~b2 

1 

b2 

. 

0 

• 

• 

0 

0 

• 

• 

a  b 

a  b  . . . 

• 

. 

b 

ml  m 

m2  m 

mm  m 

m 

For  a  basis  , the  square  array 


B 


rr 


a 


lk 


t 
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Step  II.  5:  Recover  essential  Information  of  original  tableau 


ix 


IT  +  e1 


X 


d  31) 

To  Initiate  cycle  1  the  "Information  matrix"  la  now j 


1 

0 

0  1 

1 

1  •  e  •  e 

...  1 

0 

0 

1 

0  1 

1 

°1 

c  2  » • » • 

•  •  •  c 

m 

0 

0 

0 

1 

*1 

1 

X  sees 

i 

0 

0 

0 

0  1 

■4)  j  •  •  •  • 

. . .  - b  ^ 

h 

0 

0 

0  1 

■*2 

1—  bg  . . . . 

...  — bg 

b8 

0 

e 

e 

• 

0 

1 

1 

0 1 

1 

-b 

e 

• 

• 

"fc  sees 

...  1“  b 

e 

e 

• 

b« 

m 

lu 

m 

in 

r-2o:6 
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In  inverse  form  this  becomes: 

Cycle  1: 
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4.  ip;  I  704  Pro -ram 

The  algorithm  has  been  coded  for  IBM  704.  A  brief  description 
of  the  program  and  results  will  be  given  In  this  section.  The 
code  was  written  In  the  FORTRAN  language.  It  consists  of  the 
SHARE  FORTRAN  code,RSMl,  for  linear  programming  with  3ome 
changes  and  deletions  and  the  addition  of  the  following  eight 
subroutines : 

A1  —  compute  c  P, 

A2  —  compute  solution  x  *=  ]F 

A4  —  Price  out  a  s  that  dropped  out  of  the  basis 

AXXJ  -  compute  aP  =  AxP 
BOS  —  main  routine 
CYC  -  do  one  Iteration 
NEW  —  start  phase  one 
UNL  —  unload  data  or  restart 


RSM1  performs  the  revised  simplex  method  with  the  explicit 
form  of  the  In  se.  Since  the  algorithm  for  the  chemical 
equilibrium  code  requires  the  computation  of  n  exponential 
functions  for  each  a1'  generated,  the  inverse  is  kept  in 
double  precision  form  and  the  pivoting  operation  is  done  In 
double  precision  to  preserve  accuracy. 

One  technique  used  In  3tep  I  of  the  Iterative  procedure 
for  preserving  accuracy  Is  to  compute  quantities  y'  y'  instead 

J 

of  yi;  y.  Letting  J  =  *  be  defined  by  y#  =  Max  y  ^  ,  we  set 

y\  =  y</y*»  y  =  Z  y'y  J  ~  Thus 


ir 

« 
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y*  =  1#  Yj  Is  between  0  and  1 ,  1  <(  y  n  and  In  y  =  In  y'  -  c*. 
Letting  c*  =  -in  y  # ,  then  in  y^  =  in  yj  —  c*  and  in  y  «*  in  y'  -  c# 

Instead  of  creating,  at  each  iteration,  a  new  column  a' 
ao  outlined  in  section  2,  it  is  possible  to  enter  into  the 
basi3  an  a  that  dropped  out  of  the  basis  of  some  previous 
Iteration  lc.  This  is  the  case  if 


(2h) 


G(uk)  - 


tt°  ak  <  0, 


Therefore,  each  generated  a  is  saved  and,  before  deciding  to 
generate  a  new  a  on  a  given  iteration,  (24)  13  computed  for  all 
old  a  s.  If  none  price  out  negatively,  step  I  is  used  for 
generating  a  new  a.  The  X  vector  corresponding  to  each  a  is 
also  saved  to  be  used,  if  necessary,  in  the  computation  of  the 
final  solution  in  step  III. 

m 

As  3hown  in  [5]  the  phase  one  procedure,  where  Y  x  . 

a  n  h  i 

is  minimized,  reduces  to  the  standard  simplex  method.  Thus, 
at  each  iteration  in  phase  one,  the  new  a  column  brought  into 
the  basis  is  the  original  column  a  ,  1=1,..., m,  which  has  the 
most  negative  reduced  co3t,  c..  The  X  vector  corresponding  to 

J 

tnlo  a  is  (0,  0, . . .  ,  1 ,  . . , ,0) ,  where  the  one  appears  in  column  J. 

The  program  was  written  to  handle  chemical  equilibrium 
problems  with  partitions  (chemical  phases).  The  twelve  equation 
system  representing  the  external  respiratory  system  is  such  a 
problem  (see  [6]).  The  thirty  columns  are  divided  into 


three  partitions  and  each  Iteration  Involves  generating,  from 
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one  of  the  partitions,  an  a  column  for  the  basis.  In  general, 

a  generated  a  cores  from  one  of  the  K  partitions,  where  the 

partitions  are  selected  In  cyclic  order.  The  problem  terminates 

when,  for  K  consecutive  Iterations,  =  1  in  each  partition  1. 

It  may  be  that  It  Is  not  possible,  using  the  columns  of  a 

selected  partition,  to  generate  an  a  for  the  basis  with  negative 

7.  In  ouch  a  case  the  partition  is  never  looked  at  again. 

If  each  original  column  J  13  made  a  partition,  then  x,  =  x 

In  each  partition,  and  the  logarithm  term  (In  x  /x)  drops  out 

J 

of  the  objective  function.  Thus  the  problem  becomes  a  standard 
linear  programming  one. 

The  program  allows  for  problems  up  to  dimensions  50  x  200. 
Hie  12  x  JO  problem  took  23O  Iterations  (about  minutes  on 
the  704)  to  reach  a  final  solution.  The  number  of  Iterations 
depends,  of  course,  on  the  various  tolerances  used  In  the  code. 
The  most  critical  tolerance  is  that  used  for  testing  y  =  1. 

In  the  above  example  it  was  10  . 

The  code  is  available  from  The  RAND  Corporation. 
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